*************************************************
*				sun_ab_ado.do					*
*												*
*	This do-file is called to implement the		* 
*	Sun and Abraham (2020) event-study regs		*
*	for the main and Appendix results in 		*
*	Ricks  ® Kay "Time-Limited Subsidies"		*
*												*
*												*
*												*
*************************************************

* Arguments that the do-file takes:

	* Local 1: Outcomes (main first)
	* Local 2: Heterogeneity
	* Local 3: List of Cohort-by-Event Regressors
	* Local 4: Fixed Effects
	* Local 5: Cluster


qui{

* Count the heterogeneity cuts	
if "`2'" == "" {
	local i = 0
	foreach y in `1' {
	qui di `i++'
	}

	* Set up empty matrices
	mat betas = J(`i',3,.)
	mat ses = J(`i',3,.)
	mat beta_het = J(`i',4,.)
	mat se_het = J(`i',4,.)	
	
	* Loop over outcomes
	local j = 1
	foreach outcome in `1' {

	  * Regress outcome on all cohort-by-event indicators
	  reghdfe `outcome' `3', a(`4') cluster(`5')
	
	  noi di "`outcome' regression"
	  noi di "Observations: "e(N)
	  
	  * For the "Main" Outcome calculate all ES coefficients
	  if `j'==1 {
		
		* Set up macros
		mat ES_pool = J(228,6,.)
		local big_list_short = ""
		local big_list_long = ""

		* Denominator for short-term effect
		count if count>120 & count<=144
		local big_den_short=`r(N)'	

		* Denominator for long-term effect
		count if count>144 & count<=180
		local big_den_long=`r(N)'	

		* Loop over each event period (other than 120)
		forval n = 61/180 {	
		  if `n'!=120{
			
			* Reset local 
			local list ""	

			* Denominator for each event period
			qui count if count==`n'
			local den = `r(N)'
			
			* Loop over vintages in each event (that exists)
			forval y = 2002/2021 {
				
			  * Get denominator for non-empty vinatage-event pairs	
			  cap confirm variable e`n'_c`y' 
			  if !_rc {	
				qui count if count==`n' & vintage==`y'
				
				* Vintage share
				local _p = `r(N)'/`den'

				* Define weighted average for the linear combination
				if "`list'" =="" {
					local list "`_p'*e`n'_c`y' "
				  }
				else {
					local list "`list' + `_p'*e`n'_c`y' "	
				  }
				}
				
			}
			
			* Get average and se for each event coef
			qui lincom 	`list'
			
			* Add to matrix
			mat ES_pool[`n',1] = [`n',`r(estimate)',`r(estimate)'+1.96*`r(se)',`r(estimate)'-1.96*`r(se)'] 
			}
			
			* Put zero for 120
			else{
				mat ES_pool[`n',1]=[`n',0,0,0]
			}
			
			* Calculate short-term effect
			if `n'>120 & `n'<=144 {
			  
			  * Vintage-event share within whole short-term
			  local _q = `den'/`big_den_short'
			
			  * Linear combination for short-term
			  if "`big_list_short'" ==""{
					local big_list_short "`_q'*(`list') "
				  } 
			  else {
					local big_list_short "`big_list_short' + `_q'*(`list')"	
				  }
			  }
			
			* Calculate long-term effect
			if `n'>144 & `n'<=180 {
				
			  * Vintage-event share within whole long-term	
			  local _q = `den'/`big_den_long'
			
			  * Linear combination for long-term
			  if "`big_list_long'" ==""{
					local big_list_long "`_q'*(`list') "	
				  }
			  else {
					local big_list_long "`big_list_long' + `_q'*(`list')"	
				  }
			}
			
			
		}
		
		* Get the average and se on the short-term effect
		lincom 	`big_list_short'
		
		* Get the average and se on the long-term effect
		lincom 	`big_list_long'
		}
		
	* Note: lincom can't handle the full list, so do it by hand

	
	* Reset macros

	cap mat drop P
	cap mat drop P_short
	cap mat drop P_long

	local row = 1
	
	* Consturct a vector of weights (in the post period)
	count if count >120 & count<=180
	local den = `r(N)'

	count if count >120 & count<=144
	local den_short = `r(N)'

	count if count >144 & count<=180
	local den_long = `r(N)'
		  
	* Loop over all vintage-event groups	  
	foreach v of varlist `3' {
		
		* Vintage
		local y = real(substr("`v'",strpos("`v'","c")+1,.))
		
		* Count
		local n = real(substr("`v'",strpos("`v'","e")+1,strpos("`v'","c")-3))

		* Weight is the share of observations pertaining to the main estimate
		if `n' > 120 & `n'<=180{
		  qui count if count == `n' & vintage == `y'
		  mat P = [nullmat(P) \ `r(N)'/`den'] 
		}
		
		* Weight is zero for pre period, and "min" and "max"
		else {
		  mat P = [nullmat(P) \ 0]
		}	
		
		* Weight is the share of observations pertaining to the short-term estimate
		if `n' > 120 & `n'<=144{
		  qui count if count == `n' & vintage == `y'
		  mat P_short = [nullmat(P_short) \ `r(N)'/`den_short'] 
		}
		
		* Weight is zero for pre period and long period
		else {
		  mat P_short = [nullmat(P_short) \ 0]
		}	
		
		* Weight is the share of observations pertaining to the long estimate
		if `n' > 144 & `n'<=180{
		  qui count if count == `n' & vintage == `y'
		  mat P_long = [nullmat(P_long) \ `r(N)'/`den_long'] 
		}
		
		* Weight is zero for pre period and short period
		else {
		  mat P_long = [nullmat(P_long) \ 0]
		}	
		
	}
	
	* Add one last zero for the constant
	mat P = [P \ 0]
	mat P_short = [P_short \ 0]
	mat P_long = [P_long \ 0]

	* Average Effect
	mat betas[`j',1] = e(b)*P
	mat betas[`j',2] = e(b)*P_short
	mat betas[`j',3] = e(b)*P_long

	* Get Standard Errors by Definition of Covariance
	
	* Weights
	mat W= P*P'
	mat W_short = P_short*P_short'
	mat W_long = P_long*P_long'

	* SE are product of weights and Variance-Covariance Matrix
	mata st_matrix("se", sqrt(sum(st_matrix("W"):*st_matrix("e(V)"))))
	mata st_matrix("se_short", sqrt(sum(st_matrix("W_short"):*st_matrix("e(V)"))))
	mata st_matrix("se_long", sqrt(sum(st_matrix("W_long"):*st_matrix("e(V)"))))

	* Standard Errors
	mat ses[`j',1]=[se,se_short,se_long]

	* Display Results
	noi di "Main Effect on `outcome' is " betas[`j',1] 
	noi di "With Standard Error (" ses[`j',1]

	* Drop matrices 
	foreach mat in W W_short W_long P P_short P_long {
		mat drop `mat'
	}

  ***
  * Calculate Vintage Heterogeneity	
  ***
  
    * Reset macros
	local row = 1
	cap mat drop P_2
	cap mat drop P_3
	cap mat drop P_4
	
    * Generate vintage-group weights
	count if count >120 & count<=144 & inlist(vintage,2002,2003,2004,2005,2006)
	local den_2 = `r(N)'
		  
	count if count >120 & count<=144 & inlist(vintage,2007,2008)
	local den_3 = `r(N)'

	count if count >120 & count<=144 & inlist(vintage,2009,2010,2011)
	local den_4 = `r(N)'	  
	
	* Loop over 
	foreach v of varlist `3' {
		local y = real(substr("`v'",strpos("`v'","c")+1,.))
		local n = real(substr("`v'",strpos("`v'","e")+1,strpos("`v'","c")-3))


		* Weight is the share of observations 
		if `n' > 120 & `n'<=144{
		  qui count if count == `n' & vintage == `y'& inlist(vintage,2002,2003,2005,2006)
		  mat P_2 = [nullmat(P_2) \ `r(N)'/`den_2'] 
		  
		  qui count if count == `n' & vintage == `y' & inlist(vintage,2007,2008)
		  mat P_3 = [nullmat(P_3) \ `r(N)'/`den_3'] 
		  
		  qui count if count == `n' & vintage == `y'& inlist(vintage,2009,2010,2011)
		  mat P_4 = [nullmat(P_4) \ `r(N)'/`den_4'] 
		  
		}
		
		* Weight is zero for pre period and long period
		else {
		  mat P_2 = [nullmat(P_2) \ 0]
		  mat P_3 = [nullmat(P_3) \ 0]
		  mat P_4 = [nullmat(P_4) \ 0]
		}	
		
		
	}
	
	forval num = 2(1)4{
		
		* Add one last zero for the constant
	    mat P_`num' = [nullmat(P_`num') \ 0]
		* Average Effect
		mat beta_het[`j',`num'] = e(b)*P_`num'
		* Get Standard Errors by Definition of Covariance
		mat W_`num' = P_`num'*P_`num''
		mata st_matrix("se_`num'", sqrt(sum(st_matrix("W_`num'"):*st_matrix("e(V)"))))
		
	}
	
	mat se_het[`j++',2]=[se_2,se_3,se_4]
	
	forval num = 2(1)4{
		mat drop W_`num'
		mat drop P_`num'
		mat drop se_`num'
		
	}
	
	}

	
	}
	
***	
*For Heterogeneity Casae
***

else{
	
	* Make sure splitting variable is binary
	assert inlist(`2',0,1)
	
	* Set macros
	local i = 0
	
	foreach y in `1' {
	  qui di `i++'
	  }

	mat betas = J(`i',6,.)
	mat ses = J(`i',6,.)

	local j = 1
	
	* Loop over outcomes
	foreach outcome in `1' {

	  * Regress 
	  reghdfe `outcome' `3', a(`4') cluster(`5')

	  * make  matrix
	  mat ES_pool = J(228,8,.)

	  * Loop over event periods (other than 120)
	  forval n = 61/180 {
		if `n'!=120{
			
		* Set macros	
		local list0 ""	
		local list1 ""	

		qui count if count==`n' & `2'==0
		local den0 = `r(N)'
		
		qui count if count==`n' & `2'==1
		local den1 = `r(N)'
	
		* Loop over vintages and groups to get weights
		forval y = 2002/2021 {
		  forval x = 0/1 {
		  cap confirm variable e`n'_c`y'_f`x' 
		  if !_rc {	
			qui count if count==`n' & vintage==`y' & `2'==`x'
			local p`x' = `r(N)'/`den`x''
			
			if "`list`x''" =="" {
				local list`x' "`p`x''*e`n'_c`y'_f`x' "
			  }
			  else {
				local list`x' "`list`x'' + `p`x''*e`n'_c`y'_f`x'  "	
			  }
			}
		  }
		}
		
		* Calculate event coefs for each group
		
		if 	"`list0'" != "" {
		  qui lincom 	`list0'
		  mat ES_pool[`n',1] = [`n',`r(estimate)',`r(estimate)'+1.96*`r(se)',`r(estimate)'-1.96*`r(se)'] 
		}
		
		else{
		  mat ES_pool[`n',1] = [`n',0,0,0] 
		}
		
		if 	"`list1'" != "" {
		  qui lincom 	`list1'
		  mat ES_pool[`n',5] = [`r(estimate)',`r(estimate)'+1.96*`r(se)',`r(estimate)'-1.96*`r(se)'] 
		}
		else {
		  mat ES_pool[`n',5] = [0,0,0] 
		}
		
		}
		else{
			mat ES_pool[`n',1]=[`n',0,0,0,0,0,0]
		}
		
		
	}

	* Consturct a vector of weights (in the post period)

	local row = 1

	forval x = 0/1{
	  cap mat drop P`x'
	  cap mat drop P`x'_short
	  cap mat drop P`x'_long
	
		
	  count if count >120 & count<=180 & `2'==`x'
	  local den`x' = `r(N)'

	  count if count >120 & count<=144 & `2'==`x'
	  local den`x'_short = `r(N)'

	  count if count >144 & count<=180 & `2'==`x'
	  local den`x'_long = `r(N)'
	}

	* For each vintage group event cell get the weight
	local pp = 1
	foreach v of varlist `3' {
		local x = real(substr("`v'",strpos("`v'","f")+1,.))
		local y = real(substr("`v'",strpos("`v'","c")+1,4))
		local n = real(substr("`v'",strpos("`v'","e")+1,strpos("`v'","c")-3))
		
		* Weight is the share of observations pertaining to the main estimate
		if `n' > 120 & `n'<=180  {
		  qui count if count == `n' & vintage == `y' & `2'==`x'
		  mat P`x' = [nullmat(P`x') \ `r(N)'/`den`x'' ] 
		}
		* Weight is zero for pre period, and "min""max"
		else {
		  mat P`x' = [nullmat(P`x') \ 0 ]
		}	
		
		* Weight is the share of observations pertaining to the short estimate
		if `n' > 120 & `n'<=144{
		  qui count if count == `n' & vintage == `y' & `2'==`x'
		  mat P`x'_short = [nullmat(P`x'_short) \ `r(N)'/`den`x'_short'] 
		}
		* Weight is zero for pre period and long period
		else {
		  mat P`x'_short = [nullmat(P`x'_short) \ 0]
		}	
		
		* Weight is the share of observations pertaining to the long estimate
		if `n' > 144 & `n'<=180{
		  qui count if count == `n' & vintage == `y' & `2'==`x'
		  mat P`x'_long = [nullmat(P`x'_long) \ `r(N)'/`den`x'_long'] 
		}
		* Weight is zero for pre period and short period
		else {
		  mat P`x'_long = [nullmat(P`x'_long) \ 0 ]
		}	
		
	}
	
	forval x = 0/1{
		* Add one last zero for the constant
		mat P_`x' = [P0*(1-`x')\ P1*(`x') \ 0]
		mat P_`x'_short = [P0_short*(1-`x')\ P1_short*(`x') \ 0]
		mat P_`x'_long = [P0_long*(1-`x')\ P1_long*(`x') \ 0]
		
		* Average Effect
		mat betas[`j',1+3*`x'] = e(b)*P_`x'
		mat betas[`j',2+3*`x'] = e(b)*P_`x'_short
		mat betas[`j',3+3*`x'] = e(b)*P_`x'_long

		* Get Standard Errors by Definition of Covariance
		mat W`x'= P_`x'*P_`x''
		mat W`x'_short = P_`x'_short*P_`x'_short'
		mat W`x'_long = P_`x'_long*P_`x'_long'

		mata st_matrix("se", sqrt(sum(st_matrix("W`x'"):*st_matrix("e(V)"))))
		mata st_matrix("se_short", sqrt(sum(st_matrix("W`x'_short"):*st_matrix("e(V)"))))
		mata st_matrix("se_long", sqrt(sum(st_matrix("W`x'_long"):*st_matrix("e(V)"))))

		mat ses[`j',1+3*`x']=[se,se_short,se_long]

		

	* end `x' loop
	}
	forval x = 0/1{
	foreach mat in W`x' W`x'_short W`x'_long P`x' P`x'_short P`x'_long P_`x' P_`x'_short P_`x'_long {
			mat drop `mat'
		}
	}
	
	* end outcome loop
	}
* end heterogeneity "else"
}

}

